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The baryonic acoustic signature in the large-scale clustering pattern of galaxies has been detected 
in the two-point correlation function. Its precise spatial scale has been forwarded as a rigid-rod ruler 
test for the space-time geometry, and hence as a probe for tracking the evolution of Dark Energy. 
Percent-level shifts in the measured position can bias such a test and erode its power to constrain 
cosmology. This paper addresses some of the systematic effects that might induce shifts: namely 
non-linear corrections from matter evolution, redshift space distortions and biasing. We tackle 
these questions through analytic methods and through a large battery of numerical simulations, 
with total volume of the order ~ 100 [Gpc 3 h~ 3 ]. A toy-model calculation shows that if the non- 
linear corrections simply smooth the acoustic peak, then this gives rise to an 'apparent' shifting 
to smaller scales. However if tilts in the broad band power spectrum are induced then this gives 
rise to more pernicious 'physical' shifts. Our numerical simulations show evidence of both: in real 
space and at z=0, for the dark matter we find percent level shifts; for haloes the shifts depend 
on halo mass, with larger shifts being found for the most biased samples, up to 3%. From our 
analysis we find that physical shifts are greater than ~ 0.4% at z = 0. In redshift space these 
effects are exacerbated, but at higher redshifts are alleviated. We develop an analytical model to 
understand this, based on solutions to the pair conservation equation using characteristic curves. 
When combined with modeling of pairwise velocities the model reproduces the main trends found 
in the data. The model may also help to unbias the acoustic peak. 

PACS numbers: 



I. INTRODUCTION 

Within the last few years the discipline of physical cos- 
mology has greatly benefited from a considerable influx 
of extremely high fidelity data-sets, which have enabled 
measurements of the large scale structure of the Uni- 
verse to be made with unprecedented precision; and to- 
gether these data have led to the establishment of the 
'standard model of cosmology': the flat, Dark Energy 
dominated collisionless Cold Dark Matter (CDM) model 
[1-7]. Whilst the CDM particles are well founded from 
a particle physics point of view, the Dark Energy may 
arise through a number of possible mechanisms, most 
of which are of deep consequence to much of physics if 
found to be true [8-10]. The task of modern theoretical 
and observational cosmology, therefore, is to construct 
robust tests to expose the true physical character of the 
Dark Energy and hence differentiate between hypothe- 
ses. A number of experiments are currently underway 
with this sole purpose in view, and many more are being 
planned for the future (see [9, 10] and references there 
in for a comprehensive review of current and future mis- 
sions). The Dark Energy tests fall into two main classes: 
those which perform geometric tests of gravity and those 
which perform growth of structure tests. The geometric 
tests are essentially the use of 'standard candles' (Type 
la Supernova) and 'standard rods' (baryonic acoustic os- 
cillations), whereas the growth of structure tests, exam- 
ine how the growth rate of perturbations changes as a 
function of cosmological epoch. Weak lensing by large 



scale structure and the multiplicity function of clusters 
fall into both categories and therefore potentially offer 
the most powerful discriminatory means. However, in 
order to make precise, accurate and useful constraints 
on the Dark Energy, the systematics of each experiment 
must be fully understood and controlled to sub-percent 
accuracy [9, 10] - the removal of 'unknown unknowns' is 
imperative. 

For instance, the standard candle measurement from 
Type la supernovae must address the issue of whether 
or not the ensemble of candles evolves with redshift, i.e. 
through mctalicity effects, or evolution of the underlying 
host galaxy properties as a function of redshift. More- 
over until the 'true' mechanism that drives the nova is 
understood, it may be the case that this potential sys- 
tematic can only be quantified and eliminated once the 
data are in hand. 

In this paper we shall restrict our attention to the 
second of the geometric tests, that is the standard rod 
measurement from the Baryonic Acoustic Oscillations 
(BAO). Like the standard candle test, this method also 
suffers from potential systematics; the three knowns in 
this case are: nonlinear mass evolution, non-linear bias 
and redshift space distortions (hereafter, we shall refer 
to these together as clustering systematics). However, 
unlike the case for Type la Supernova, because the pro- 
cesses driving any possible evolution are plausibly under- 
standable ab-initio, there is not much room for unknown 
unknowns and there is some hope for estimating and mit- 
igating these effects well-before the data streams in from 
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the next generation surveys. This is important because 
if the BAO peak is displaced by even 1%, this will induce 
a bias in the inferred value of the dark energy parameter 
w on the order of 5% [7, 11]. 

The physical picture for the BAO signature is as fol- 
lows: before the epoch of recombination, acoustic oscilla- 
tions were able to propagate through the photon-baryon 
plasma at the sound speed, and these waves were weakly 
coupled to dark matter through gravity. After recom- 
bination the photons free stream out of the perturba- 
tions and this gives rise to the observed CMB ([12]), the 
dark matter and segregated baryons then relax together 
over time and the self-same acoustic features that are im- 
printed in the CMB become imprinted in the dark mat- 
ter distribution. The characteristic scale for the acoustic 
waves is set by the sound horizon at last scattering 
(see [13] for a description of how to calculate this), and 
this in turn imprints a characteristic scale in the pattern 
of galaxies and it is supposed that this has the properties 
of a 'standard rod'. 

The BAO features have been detected by various 
groups: in the two-point correlation function of Lumi- 
nous Red Galaxies (LRG) by [7], and in the power spec- 
trum of galaxies and LRGs by [2, 6, 14-20]. The BAOs 
have also been the subject of much vigorous theoretical 
and numerical research [11, 21-30, 32-38]. The ques- 
tion of whether there are non-linear effects at play on 
the acoustic scale, is not an open question [21], however, 
whether these non-linearities give rise to an actual mo- 
tion of the acoustic peak - apparent or physical - is of 
great debate, and the most recent literature concerned 
with this question reaches conflicting conclusions: [39] 
used the fitting formula for the power spectrum from [40] 
to conclude that, there is a shift due to nonlinear mass 
evolution on the order of ~ 2% at z = 0. [23] used numer- 
ical simulations to show that there were changes to the 
broad band power spectra of dark matter and haloes, and 
in both real and redshift space, however they argued that 
provided these were accounted for, no overall shift in the 
acoustic peak position would be induced. [35] used nu- 
merical simulations with improved resolution to convinc- 
ingly confirm the results from [23], that the power spec- 
tra were not immune to strong broad band tilts. Based 
on these results they suggested that percent level shifts 
in the position of the acoustic peak were highly plausi- 
ble. The main findings of these works were most recently 
substantiated by [11]. On the other hand, [30] used a 
model based on Lagrangian displacements of the initial 
density distribution to argue that any acoustic peak shift 
in the dark matter should be only of the order 10~ 4 at 
z = 0, although they do note that "galaxy bias could 
produce a sub-percent shift" . In addition, [36] studied 
how a relatively (by BAO standards) large peak in the 
initial power spectrum evolved in numerical simulations 
and concluded that there were no noticeable shifts, in 
agreement with [30]. 

In what follows, we examine this issue in detail. We do 
this in a two-fold way: Firstly, we generate a large ensem- 



ble of large volume numerical simulations to quantify the 
possible effects. Secondly, we develop a new analytical 
model, which is based on a new solution for the pairwise 
conservation of particle pairs. When combined with a 
careful modeling of the divergence of pairwise velocities 
beyond linear theory this method is shown to capture the 
main effects that are found in the the simulations. 

The subject of this paper is therefore to answer the fol- 
lowing important questions: Does non- linear evolution 
generate a displacement of the peak of the correlation 
function? If so does the observed shift depend on the 
halo/galaxy sample considered and how? Recently, there 
has been a number of different approaches to estimat- 
ing the sound horizon from observational data, however, 
so far as we know, it has not been shown that any of 
these estimators are consistent, unbiased, or indeed min- 
imum variance estimators. The results presented in this 
work and from our previous study of the power spectrum 
[35, 37] should act as an important empirical guide for 
constructing such quantities. 

The paper is structured as follows: In Section II we 
discuss a toy model that shows that an effective smooth- 
ing of the acoustic peak in the two-point function leads 
to an 'apparent' motion of the peak. Here we also show 
how if nonlinear evolution induces a broad band tilt in 
the underlying linear power spectrum, further shifts in 
the peak position are to be expected - these we shall 
class as 'physical' shifts. Then in Section III we describe 
our ensemble of numerical simulations and present our 
measurements for the two-point correlation function of 
dark matter and haloes in real and redshift space, in- 
cluding a detailed analysis of our data. In Section IV we 
describe our new physical model and demonstrate how it 
gives rise to a transformation of the structure of the peak 
in the dark matter and halo correlation functions - and 
that this gives rise to a physical motion of the peak. We 
also compare our analytic model with the results from 
the numerical simulation and show that they are in close 
agreement. Finally, Section V summarizes our results, 
and discusses them in the wider context. 



II. APPARENT AND PHYSICAL SHIFTS 

A. Motivation 

Motivated by the calculation of the real space dark 
matter correlation function in renormalized perturbation 
theory (hereafter, RPT) [41, 42], we can write the ob- 
served correlation function in terms of the linear one 
through the following relation: 

UW « f &n(r - r') K(r') d 3 r> + £ mc (r) , (1) 

where the first term on the left-hand-side represents the 
linear correlation function convolved with some symmet- 
ric kernel, K(r), and the second term, £ mc , describes 
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any effects due to non-linear mode-coupling. The dis- 
tinction between these two terms may be more clearly 
seen in Fourier space: the first term is directly propor- 
tional to the linear power spectrum at the same scale, 
and the second term represents a weighted sum over the 
information from different neighboring wavemodes. Note 
that such decomposition can always be made. In RPT, 
the kernel K is well approximated by a Gaussian [42], a 
result that becomes exact in the Zel'dovich approxima- 
tion [30, 31, 42]. 

Setting aside £ mc for the moment, we remark that it 
is sometimes thought that convolution with a Gaussian 
does not lead to a shift in the BAO peak position. In the 
following sub-section we will show explicitly that this is 
not correct and that the convolution with a symmetric 
filter does shift the peak, and that this is solely due to the 
fact that £ii n is not symmetric about the acoustic peak. 
However, as we mention in the following sub-section this 
apparent shifting of the peak may be corrected for. 

Returning now to the issue of mode coupling, as we will 
show in this work through our numerical simulations and 
through our theoretical analysis, the term £ mc in Eq. (1) 
gives rise to an actual 'physical' shift towards smaller 
scales as the clustering evolves. For reasons which are 
now clear, we shall now refer to the shifts that are gener- 
ated by the first term as being apparent, and those due 
to the second, as being physical. In the next subsection 
we present a toy-model to further illustrate the meaning 
of these terms. 



B. A toy model for the shifts 



from T 2 (k) (this is a toy- model). Restricting the range 
of interest to be small enough so that Smooth (r) is close 
to a power-law, then we would have something like our 
Eq. (2). 

The presence of the power-law means that the loca- 
tion of the local maximum, say r max , will differ from r p . 
Requiring d^/dr = means 
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If (r p — r m ) = e <C a (meaning the offset from r p 
is small compared to the width of the bump), then this 
becomes 
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The fact that we call it a bump means that Aq > A p . In 
addition, we are interested in the case where a <c p , thus 
our final expression for the peak in the linear correlation 
function is 



Part of the following analysis was inspired by ideas 
first presented by [39]. In that work one of the issues 
addressed was the apparent shift of the acoustic peak po- 
sition, induced by an inhomogeneous selection function. 
Here we use similar arguments, but directly connected 
to the distortions induced by the non-linear clustering 
transformation and bias, to examine the apparent shifts. 
Those familiar with the analysis of [39] may wish to jump 
directly to Eq. (7), which should be familiar. 

To begin our toy-model, let us suppose that the linear 
theory correlation function can be well approximated by 
a power-law plus a Gaussian bump with peak position 



located at r p : 



£(r) = A p (J±Y + A G exp 



(r- 



2a l 



(2) 



This is a reasonable approximation, since the transfer 
function can be decomposed into a smooth component, 
which models the suppression of dark matter fluctua- 
tions due to radiation dominated growth and baryon 
drag effects, and an oscillatory piece that comes from 
the baryons clumped around the sound horizon: i.e. 
T(k) = T smooth (fc) + T B Ao(k) (sec [13]); on squaring and 
Fourier transforming we get £(r) = £ smoo th (r) + £bao(»"), 
where we have for simplicity neglected the cross-terms 
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a 



A p 
~Ag 



(7) 



This shows that the fractional shift from r p is large if 7 
is large (meaning the amplitude of the power law com- 
ponent is changing rapidly), or if a jr p is large (meaning 
the bump is broad, so the change in the amplitude of 
the power law component matters), or if A p /Ag is large 
(meaning that the power-law component matters). 

What concerns us now is: How does the peak scale 
change if one of our clustering systematics alters one or 
all of these terms? Suppose Aq — ► Aq(1 + Sa g ), A p — ► 
A p (l + 5a p ), etc., then we would have 



eiii 



(l + ^ 7 )(l + (5 CT ) 2 (l + J Ap ) 
(1 + <W) 



(8) 



and if these changes are all small, then e — ► eii n (l + S e ) 

<5 £ w 6 1 + 2S a + 5 Ap - S Ag . (9) 

If the only effect of the clustering systematics is to 
smooth out the spike to a bump, then they may simulta- 
neously increase the width of the peak and decrease A G : 
i.e. Aq oc I/it, implying that S e However, because 

5 a can be larger than ~ 0.1, the effect on e oc (1 + 5 a ) 3 
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FIG. 1: Halo correlation functions at z = in real (top) and redshift (bottom) space. Different symbols in each panel show 
results for massive (top) to less massive halos (bottom). Table I gives the precise bins in halo mass. Error bars come from 
the dispersion between the measured £ in our 50 simulations; a total volume of 105 (/i _1 Gpc) 3 . Solid line in top panel shows 
the linear theory correlation function multiplied by an arbitrary constant so that it approximately matches the signal from the 
intermediate mass bin. Vertical dashed line shows the position of the acoustic peak in this linear correlation function: it lies at 
106 h' 1 Mpc. 



may be substantial. We emphasize that such an appar- 
ent shift would occur even if there were no physical shift 
in the position of the peak. Turning now to the physical 
shifts: if J 7 7^ then we shall say that our clustering 
systematics have changed the underlying power-law and 
that this will lead to a physical motion of the acoustic 
peak. 



Before we move on, we note that there are circum- 
stances under which the apparent shifts may be consid- 
ered as benign and so removed, namely the Gaussian 
smoothing case. However, the physical shifts are more 
pernicious and when these distortions are present it is 
not clear how best to reconstruct the unperturbed peak 
for both of the shifts. We shall reserve further discussion 
of this matter for our future work. However, we note 
that for dark matter in real space these effects can be 
calculated rather precisely [37]; in particular, the physi- 
cal shifts are more complicated than just an overall tilt of 
the underlying power-law as described by our toy model. 



III. APPARENT AND PHYSICAL SHIFTS 
FROM NUMERICAL SIMULATIONS 

A. The ensemble of simulations 

For the range of cosmologies that are acceptable, the 
BAO peak is located at about r p ~ 100 h~ x Mpc. A 
large simulation volume is therefore required in order 
to minimize the cosmic variance in the measurement on 
these scales and also to correctly account for the mode- 
coupling from scales beyond r p that may drive evolution 
[35]. However, to control the sample variance down to a 
level of a few percent requires the generation of a huge 
computational volume. To make this task feasible, given 
our finite computer resources, we decided to run a large 
ensemble of large simulations as opposed to one single 
extremely large simulation. As we will show this allowed 
us to robustly answer the question as to whether there 
is any apparent or physical evolution in the peak posi- 
tion. These simulations will also allow us to assess how 
sensitive future surveys will be to measuring the acoustic 
feature. To this end, we have run fifty realizations of cu- 
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FIG. 2: Same as Fig. 1 but at z = 0.5. 



bic boxes with side Lbox = 1280ft. _1 Mpc, giving a total 
comoving volume of about 105 (h" 1 Gpc) 3 , just under 
four times the volume of the Hubble volume simulation. 
This is approximately the volume ADEPT plans to sur- 
vey, and is more than an order of magnitude larger than 
any current or proposed ground based experiment [43] . 

The cosmological parameters for the ensemble were se- 
lected to be in broad agreement with the WMAP best fit 
model [1]: Q m = 0.27, fl A = 0.73, fl b = 0.046, h = 0.72 
and as(z = 0) = 0.9. For this cosmology, linear the- 
ory predicts the position of the acoustic peak, i.e., the 
local maximum of the auto-correlation function of dark 
matter, to occur at 106 hr 1 Mpc. 

Each simulation was then run with 640 3 particles. We 
used the cmbf ast [44] code to generate the linear theory 
transfer function, and we adopted the standard parame- 
ter choices, but took the transfer function output rcdshift 
to be at z = 49. The initial conditions for each simula- 
tion were then laid down at z = 49 using the 2LPT code 
described in [45, 46] and subsequent gravitational evolu- 
tion of the equations of motion was performed using the 
Gadget2 code [47]. Each realization ran to completion in 
roughly 1900 timesteps from redshift z = 49 to z = 0, 
and the comoving force softening was set at 70 kpc/h. 

Haloes were identified in the redshift z = 0, 0.5 and 
1 outputs of each realization, using the friends-of-friends 
algorithm with linking-length parameter I = 0.2 (this 



TABLE I: Halo samples as a function of redshift. Halos in 
the "large", "intermediate" and "small" mass bins M > M3, 
M 2 < M < M z and Mi < M < M 2 , respectively. Masses 
are in units of /i _1 Mq and comoving number densities fin in 
(ft^ 1 Mpc)" 3 . 





2 = 


z = 0.5 


2 = 1 


n H 


M 3 


1.5 x 10 14 


10 i4 


5.7 x 10 i3 


1.9 x 10 _b 


M 2 


7 x 10 13 


5 x 10 13 


3.1 x 10 13 


3.4 x 10~ 5 


Mi 


4 x 10 13 


3 x 10 13 


2 x 10 13 


4.8 x 10~ 5 



choice is standard) . Halo masses were then corrected for 
the error introduced by discretization of the halo density 
structure [48] . Since the error in the estimate of the halo 
mass diverges as the number of particles sampling the 
density field decreases, we only study haloes containing 
33 particles or more. At each redshift we present results 
for the three bins in halo mass. These bins were chosen 
by counting down in mass from the most massive halo, so 
that the number in each bin is the same at each redshift. 
Table I shows the resulting cuts in halo mass, and the 
associated comoving number densities. 
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FIG. 3: Same as Fig. 1 but at z = 1. 



B. The measured correlation functions 

The correlation functions were estimated using the 
standard estimator: £(r) = DD(r) / RR(r) — 1, where 

£ is the bin averaged correlation function, DD(r) is the 
number of true data pairs in the bin and RR(r) is the 
number of pairs expected after we randomize the posi- 
tions. We also note that when dealing with the rcdshift 
space data, we apply the distortion separately in the x—, 
y— and 2— directions and measure three correlation func- 
tions, these are then averaged together to construct a 
single estimate for each realization. 

Figure 1 shows the auto-correlation functions of the 
halos in each of the selected mass bins in our 2 = out- 
puts. Top and bottom panels show £(r) and £ s (r), the 
real and redshift space correlation functions. We have 
chosen to show £(r) rather than r 2 £(r) because, as dis- 
cussed earlier, the peak in the former is more directly 
related to the sound horizon scale r s . The error bars on 
the data points come from the scatter around the mean 
value of £ as measured in the fifty realizations (i.e. from 
the diagonal elements of the covariance matrix divided 
by the square-root of the number of realizations, which 
for our case is: \/50 ~ 7). 

The solid line in the top panel shows the dark mat- 
ter correlation function predicted by linear theory, mul- 



tiplied by a constant factor so that the curve approxi- 
mately matches the signal seen in the intermediate mass 
bin on scales r < 80 ft. -1 Mpc. The vertical dashed line 
shows the location of the local maximum in this func- 
tion: 106 hr 1 Mpc. Considering the results in real space 
(top panel), the figure clearly shows that the local max- 
ima of the measured correlation functions are systemat- 
ically shifted to smaller scales compared to this mark. 
Moreover, it appears that the magnitude of the shift 
steadily increases with halo mass. Turning to the results 
in redshift space, we see that this effect is even more 
pronounced. 

Figures 2 and 3 show results at redshifts, 2 = 0.5 and 
1. Although the distortions from the linear case appear 
to be slightly smaller, we again see clearly that the trends 
are similar to those of the redshift zero case. 

We now draw attention to another point of interest. 
As is expected, these selected halo samples are signif- 
icantly more clustered than the mass. The large-scale 
bias factors, as measured by the (square root of the) ratio 
of the halo correlation function to that of the measured 
dark matter on scales ~ 70 ft. -1 Mpc (where nonlinear ef- 
fects appear to be small) are b — 1.4, 1.8, 2.6 for the 
z = halos, b = 1.9, 2.3, 3.2 for the z = 0.5 halos, and 
b = 2.5, 3.0, 3.9 for the 2 = 1 halos - with the most 
massive halos having the largest bias parameters. What 
is not so obvious now is that the halo clustering signal for 
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each bin at the three different redshifts is almost constant 
in time. For reference, consider the linear theory growth 
factor which is smaller by a factor of 0.785 between red- 
shifts z = and 0.5, so the amplitude of drops be- 
tween at z — and 0.5 by a factor of ~0.615. This result 
is a direct consequence of studying the signal at fixed co- 
moving number density: whilst the clustering of the mass 
is much smaller at higher redshift, the high redshift halos 
are significantly more biased. At fixed number density, 
the two effects approximately cancel out, keeping the net 
clustering signal fixed. This is important in view of the 
fact that galaxies of approximately constant comoving 
density represent a popular choice for the target sample 
galaxy to measure the BAO signature over a range of 
redshifts, i.e. the Luminous Red Galaxies (LRG). 



C. Statistical properties and model fitting 

In this section we examine the statistical properties of 
the halo-halo correlation functions and present our model 
fitting analysis. 

Figure 4 shows again the simulation results in real 
space and for the smallest (top) and largest (bottom) 
bins in halo mass, at z = (left) and z — 1 (right). Fig- 
ure 5 shows again the results in redshift space. In each 
panel, symbols show the mean value of £ hh for the given 
bin in r, averaged over the 50 simulations; shaded regions 
show the standard deviation over the 50 realizations, and 
error bars show the error on the mean (they are smaller 
than the shaded regions by a factor of \/50 ~ 7) . 

The first point to note is that the scatter amongst real- 
izations is remarkable, given that each one of our boxes is 
about three times larger than the volume probed by the 
SDSS LRG sample. We further emphasize that at least 
in redshift space this is likely to be a lower bound on the 
true underlying scatter, since there is no contribution 
to the variance from virial motions of the dark matter 
particles or galaxies. Clearly, enormous volumes will be 
required to measure £ hh , and thus the galaxy correlation 
function, to the precision required for percent precision 
cosmology, and this justifies our earlier assertion at the 
beginning of this Section. 

Comparing now the scatter exhibited in the z = real 
space low mass halo sample with that found for the high 
mass sample, at a first glance we see that the scatter ap- 
pears to decrease as halo mass increases; and this trend is 
also exhibited in redshift space data. However as one goes 
to higher redshifts no obvious trend is apparent between 
low and high mass samples. Comparing halo samples at 
the same fixed number density but at different epochs, we 
see that the scatter is much reduced for the low mass halo 
sample, but roughly constant for the higher mass sample. 
This suggests that what is meant by 'high' or 'low' mass 
is a very subjective quantity: 'low' mass here must mean 
relative to the typical halo mass at that epoch. However 
as we shall show shortly these trends with halo mass and 
measured epoch can not be characterized so naively 



We now turn to our modeling of the data. Based on our 
discussion from Section II, we now attempt to fit the cor- 
relation functions by assuming that each can be described 
as a linearly biased version of the linear theory correla- 
tion function, smoothed with a Gaussian filter. There 
are thus two free parameters for such fits-the scale of the 
Gaussian smoothing filter Rq and the overall amplitude 
of the bias, which we define b\ = £, hh (r)/£,Lin(r), where 
£ hh (r) is the halo-halo correlation function. Note that 
for this theoretical case we shall assume that the trans- 
fer function and hence the cosmological model are fully 
specified, which for our simulations they are of course, 
however in reality one should consider fitting for these 
parameters jointly with other cosmological parameters, 
albeit with constraints from external data i.e. CMB etc., 
since these should not be considered known a priori for 
a realistic case. 

The best-fit model parameters for each sample were de- 
termined by generating a 2D cubical grid of models over 
the ranges b x e [0.5,10.0] and R f e [0.5,10.0], spaced 
by steps of 0.01, and then computing the x 2 for each, 
with the final best-fit model being identified as the one 
with the minimum returned value. Explicitly the \ 2 wc 
minimize is 



(10) 



where we define the ensemble average difference vector 
Y = y-y moA {^\b u R f ) , (11) 

with y = (£i h , . . . , and x = (n, . . . , rjy), and where 

y mod is a vector of model values. C^hh} is our estimate of 
the covariancc matrix of the mean correlation functions 



<c hh > 



y T ®y-y T • 



y 



IN, 



real 



(12) 



where ® is the direct product operator and we divide by 
the number of realizations, l/iV rca i, because we are esti- 
mating the covariances of mean quantities. The inversion 
of the above Covariance matrix was performed using the 
SVD algorithm [55]. The models were fit over the range 
of scales (65 /i -1 Mpc < r < 140/i _1 Mpc) in 21 equal 
bins in radius giving ~ 3.5 hr x Mpc per bin. 

Consider now the three thin solid curves, the central 
line represents the best-fit smoothed and scaled linear 
theory correlation function; and the best-fit values for b\ 
and Rf, along with their respective \ 2 values are also 
reported in each panel. Interestingly, we note that the 
bias factors in the redshift space figures, recovered from 
this fitting procedure, roughly agree with our earlier es- 
timates from simply considering the data points around 
r <~ 70Mpc//i. For comparison the dashed lines show 
the best-fit unsmoothed linear theory correlation func- 
tion with a linear bias. Clearly this model does not pro- 
vide a reasonable fit to the simulation data over the range 
of scales presented. 

We may now make a prediction for the variance of 
the correlation function in each bin by assuming that 
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FIG. 4: Mean (solid points), scatter (shaded region) and error on the mean (error bars) for the halo-halo correlation functions 
measured in the ensemble of 50 simulations. The long dashed curves show the linearly biased, linear theory; the central solid 
curve shows linear theory, smoothed with a Gaussian filter radius R and linearly biased b (best fit values for these parameters 
are expressed in the figure annotations). The inner and outer solid curves enclosing the best fit model show the expected 
scatter in the continuum limit and the discrete Poisson sample limit, respectively - see text for full explanation. The vertical 
lines represent the local maximum of the linear theory £ (right most dash line) and the best fit smoothed linear theory model 
(solid line) and the best-fit Tchebychev polynomial fit (triple dot dashed lines). The bottom panels show the ratio of the 
measurements to the central solid line and again the error bars are the errors on the mean. 
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FIG. 5: Same as previous figure, but in redshift space. 
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the underlying density field from which the correlation 
functions were generated is well described by a Gaussian 
random field. Following [49-53], the covariance matrix 
for bin averaged correlation function, in the limit of small 
bin sizes Ar/r <C 1, can be written 



<cer>-<e>c> 

1 f dkk 2 

(k) ,(13) 



contributions to the variance from the connected trispec- 
trum and bispectrum, [52-54] and that the discrepancy 
between high and low-mass halo samples may be caused 
by the effect of non-linear halo bias terms entering the 
variance estimate. 

Returning to the issue of how the scatter depends on 
the samples. Considering again our expression for the 
Gaussian error on the power spectrum, Eq. (14), we see 
that the relative error per mode can be written: 



where £f h is the bin averaged correlation function in bin 
i, Vfj, is the simulation volume, jo = sin(a:)/a; is the ze- 
roth order spherical Bessel function and where crp hh is 
the Gaussian variance per mode in the halo-halo power 
spectrum, which for a discrete Poisson sampling of the 
halo field, can be written 



hh 

phh v A 



n H b 2 Pun(k,z) 



(17) 



2 

O phh 



(k) 



P hh (k) + — 

n H 



(14) 



where we write our smoothed halo-halo power spectrum, 
at linear order in the over-density perturbation and bias, 
as 

P hh (k\R f ) = V,(\5 h (k\R f )\ 2 ) = blP Un (k)W 2 (kR f ) . 

(15) 

For the purposes of numerical evaluation of the above 
formulae, we follow [50] and note that one may rewrite 
the contribution to the covariance coming from the term 
l/n% as follows: 



Thus we see that the relative scatter may increase in three 
ways: as halo bias decreases; as halo number density de- 
creases; and as the power spectrum amplitude changes 
with time. For our constructed samples these effects con- 
spire in such a way that it is no longer trivial to isolate 
trends. Rather than attempting to follow this path we 
simply note that to say anything substantiative we must 
consider the full covariance matrix, since a decrease in 
diagonal covariance can be traded for an increase in off- 
diagonal covariance - which is just as important in the 
fitting. We shall reserve the important issue of power 
spectrum and correlation function covariance for a future 
study. 

Lastly, we note that for very sparse samples of haloes 
the theoretical prediction given by Eq. (13) actually 
represents a lower bound, since there should be a fur- 
ther (non-Gaussian) shot-noise contribution of the order 
S^2^ h / [njiVftV^i)] [50]. However, for our purposes 



2 , d 3 k 2 "ij^i I Vh "m v v i)j L"»J- 

— / 77TYtjo{kri) jo(krj) = <$f _ 2 , (16) this term is unimportant, since ^ h < 1. 



V, 



where the volume associated with each shell is: V{r{) = 
AirrfAr. The variance in the correlation function is then 
simply given by setting i = j in Eq. (13). Using our 
Gaussian smoothed linear model for P hh (k) ensures that 
the integrals over the Bessel functions converge rapidly. 
For a more detailed discussion of convergence properties 
sec [50]. 

Considering again the best-fit smoothed model in each 
panel, surrounding it are two sets of solid lines, a thick 
inner set and a thin outer set. The inner lines show the 
scatter between realizations that one would predict using 
the continuum limit of Eq. (14), that is when 1/uh — * 0. 
In this case, the theoretical predictions clearly underes- 
timate the true scatter in £ hh for all bins in halo mass, 
with the discrepancy being slightly worse for the lowest 
mass bin. The outer solid curves now show the effect 
of including the discreteness contribution from 1/nn in 
Eq. (14). Note that in implementing Eq. (16) it was es- 
sential to correctly account for the binning in the correla- 
tion function. In most cases, and especially for the z — 1 
correlations, this additional contribution provides a much 
improved description of the measured scatter. However, 
at low redshift and for the lowest mass haloes considered 
the theoretical estimate of the scatter appears too small. 
This is likely a consequence of the growing non-Gaussian 



D. Evidence for shifts 

We now illustrate very clearly the effects of 'apparent' 
shifts on the peak position of the correlation function, 
arising from the operation of smoothing, and explore the 
hypothesis that the peak position also exhibits large-scale 
'physical' motion. 

Before we commence with this investigation we note 
that since our simulations used the transfer function gen- 
erated at z = 49, our linear theory acoustic peak is 
slightly less sharp than that which obtains from using the 
transfer function at z = (which is more correct). This 
means that the apparent shifts will be overestimated, 
since, as was discussed in Section II, smoothing affects 
more a weaker peak. For dark matter clustering, a calcu- 
lation using renormalized perturbation theory [37] shows 
that the apparent shift is overestimated by about a fac- 
tor of two by using the transfer function calculated at 
z = 49. However, the physical shifts are not changed 
significantly. 

Considering again Figs 4 and 5, as noted above, the 
(blue) dashed lines in each panel show the associated bi- 
ased but unsmoothed linear theory correlation functions, 
and clearly these do not provide a good fit to the data. 
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TABLE II: Shifts in the BAO peak position as a function of halo sample in real space. 

SAMPLE Rfjh^Mpc] h 8 Sp [%] S Tp [%] 5 Phys [%] ~ 

2 = 2 = 1 2 = 2=1 Z = 2=1 2 = 2 = 1 2 = 2 = 1 

Mi 5.35 2.52 1.41 2.32 1.5 0.00 1.10 0.85 0.38 0.85 

M 2 5.25 2.52 1.71 2.72 1.32 0.00 5.20 0.09 2.92 0.09 

M 3 6.51 3.94 2.52 3.94 3.30 0.47 4.40 0.57 1.13 0.09 



TABLE III: Shifts in the BAO peak position as a function of halo sample in redshift space. 

SAMPLE Rflh^Mpc] b 1 S Sp [%] S Tp [%] S Phys [%] 

2 = 2 = 1 2 = 2 = 1 2 = 2 = 1 2 = 2 = 1 2 = 2 =1 

Mi 6.26 3.38 1.56 2.57 2.83 0.18 4.90 2.16 2.07 1.98 

M 2 6.97 4.04 1.91 2.88 4.60 0.66 6.13 1.20 1.51 0.57 

M 3 6.92 4.54 2.62 4.09 4.40 0.70 NA 1.80 NA 1.10 



The corresponding vertical dash lines show the position 
of the local maximum of this curve - the unperturbed 
acoustic peak: r\j p = 106 h^ 1 h^ 1 Mpc. Now consider 
the best fit smoothed models (central solid magenta), 
these clearly (by eye) provide a much improved fit to the 
data. On finding the local maximum of these smoothed 
models, we see that in all cases the peak has been shifted 
to smaller scales. These are denoted in each panel by 
the solid magenta vertical lines and with their values re- 
ported in the top right of each panel as r Sp , with 'Sp' 
meaning smoothed peak. 

Considering these smoothed model inferred peak posi- 
tions, we note several important trends: 

• All maxima lie on smaller scales than nj p ; 

• For halo samples considered at the same epoch and 
in both real and redshift space, the shifts from r\j p 
increase with increasing halo mass; 

• Considering halo samples of the same fixed number 
density at different rcdshifts, the shifts are reduced 
for the higher redshift samples; 

• Shifts are increased in the redshift space; 

• Best fit filter scale increases with halo mass. 

The spread of values in the smoothed model shifts 6s p = 
[nj P — ^SplAup, span the range 5s P ~ 1.0% — 5%, with 
the largest values being obtained as per the trends de- 
scribed above. These shifts, it can be argued [30], fall 
under the banner of apparent shifts - arising from the 
local collapse and rearrangement of matter. However, 
following our discussion of the transfer function, we ex- 
pect that these shifts would be reduced by a factor of 
~ 2 for the z = transfer function: <5s p ~ 0.5% — 2.5%. 
We also note that in the recent literature a number of 
procedures have been proposed to tackle these apparent 
shifts and, modulo the choice of filter function, most of 
these methods should be able to successfully correct for 



these effects. However, we now draw attention to the last 
of our bullet points and note the fact that the best fit fil- 
ter scale Rf, increased with halo mass. This implies that 
methods that are tested and tuned to extract BAO infor- 
mation using only the dark matter distribution will fail 
to incorporate this effect - we return to this in Sec. Ill E. 

Turning now to the question of 'physical' shifts, it is 
clear that the smoothed model does not provide consis- 
tently good fits to the measurements for all our samples. 
To see this more clearly, the bottom section of each panel 
shows the ratio of the measured points to the best-fit 
smoothed linear model. From examination of these re- 
sults it is clear that there is some evidence for structure 
in these residuals — typically, on scales smaller than the 
true acoustic peak position we find that the data points 
lie above the best fit smoothed model. This trend is 
most apparent for the present day high mass halo sam- 
ples, but is less clear for the lowest mass. This can be 
further quantified by use of the x 2 test as an indicator for 
the 'goodness-of-fit': for N = 21 data bins and m = 2 pa- 
rameters, the probability P(x 2 > 36.19|n = 19) = 0.01. 
Thus on inspection of the \ 2 values in Figs 4 and 5 we 
see that only in two instances are the mean data in agree- 
ment with this and these are for the present day low mass 
samples in real and redshift space. Based on these data 
we are led to reject our null hypothesis and accept the 
possibility that there is a physical motion of the peak. 

One alternative to the 'physical motion' hypothesis is 
that the filter choice is somehow special - and had we 
chosen the 'special' filter then this would reconcile our 
results. This view is problematic, since in order to model 
all of the above trends such a filter would have to be 
highly contrived. Thus, based on the above evidence, it 
seems that something like the second term in Eq. (1) is 
present and generating a shift in the position of the peak. 

In the next section, we provide details of a physically 
motivated model that may offer some insights into the 
origin of the dependence on halo mass of the shifts in the 
acoustic peak position. 
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Before continuing, it is of interest to characterize the 
peaks in the correlation function data using a purely arti- 
ficial parameterized model that simply matches the data 
in the least square sense. For this we write the model 
of the correlation function as a sum over the orthogonal 
Tchebyshev polynomials, Ti(x), i.e. 

m 

2/mod(r) = £ hh (r) = <HTi[x(r)] ; (18) 

i=0 

and x{r) here is a mapping that transforms the r-range 
into the range x G [—1, 1] so that we may use the normal- 
ized polynomials. We then, as before, find the coefficients 
di that minimize our % 2 Eq. (10) and polynomials up to 
degree 9 were used to describe the data and, for simplic- 
ity, in the fitting we have used only the diagonal elements 
of the covariance matrix to make the constraints. These 
artificial models along with the locations of their local 
maxima at the acoustic scale are presented in Figs 4 and 
5 as the (green) triple-dot dash curves, and the values 
of the maxima are noted in the top-right hand corner of 
each panel as rx p . For the high mass samples the shifts, 
<$t p = [njp - ?"Tp]/njp, are significantly enhanced relative 
to linear model, and for the case of the highest mass halo 
sample at z = in redshift space, the best fit model had 
no local maximum. 

We may now be more definite in what we mean by a 
physical shift: we shall say the percentage physical shift 
away from the true peak is <5ph ys = \$Sp — <>Tp|- The phys- 
ical shifts appear to be smaller than the apparent shifts 
and are roughly of the order ~ 0.4 — 3.0% at z = 0.0 
for real and redshift space data and they are somewhat 
mitigated at higher redshifts. Tables II and III collect 
together the apparent, total and physical shift values for 
the halo samples in real and redshift space, respectively. 
As we shall discuss in the following section these physical 
shifts are not accounted for in the recent methods pro- 
posed to correct for the non-linear evolution of the BAO 
peak position. 

E. Alternative BAO reconstruction methods 

Recently a number of procedures have been proposed 
to correct the BAO peak position for the effects of 
the large scale structure systematics (see for example 
[11, 16, 23, 28, 30]). These methods essentially allow 
for smoothing in configuration space and possible broad 
band tilts in the underlying power spectrum. Since the 
correlation function is the Fourier space dual of the power 
spectrum, these methods should also be equally applica- 
ble in configuration space. 

Firstly, consider those methods that attempt to correct 
the measured power spectrum for tilting by fitting away 
an arbitrary constant - this is in response to ideas from 
the Halo Model that lead us to consider a generalized 
Poisson shot-noise correction (see for example [56]). On 
Fourier transforming this model it does nothing to the 
peak of the correlation function and so may be dismissed. 
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FIG. 6: Expected behavior of the residuals in the correla- 
tion function arising from the alternative approach to BAO 
model fitting of [30], where the broad-band, but smooth, lin- 
ear power is restored. From top to bottom the lines show the 
effects where a smoothing scale of Rf — {2, 4, 8} h~ 1 Mpc/h 
were taken. This method does not replicate the S shaped 
structure of the residuals found in the data. 



Secondly, consider a model that is designed to damp 
out the acoustic oscillations, but then restore the linear 
theory power using a smooth (no BAO) version of the 
linear power [30], e.g. 

P NL (k\R) = e- k ' 2R2 b 2 Pun(k) + (l-e~ k2R2 )b 2 Pl™ oth (k), 

In the configuration space the first term transforms into 
the smoothed linear model Eq. (15), and as we saw 
this will generate apparent shifts. Considering the sec- 
ond term, we see that this function has no information 
about the acoustic scale. Fig. 6 shows that this model 
for a range of smoothing filter scales - in all cases, it 
is flat around the acoustic peak. Moreover, the ratio 
^Nh{r\R)/^un(r\R) < 1 for all smoothing scales, whereas 
the measured residuals can exceed unity on scales smaller 
than the acoustic peak (cf. Figures 4 and 5). We there- 
fore deduce that a model like equation (19) is inadequate. 

Lastly, we note that a more sophisticated method for 
correcting the spectrum for the non-linear systematics 
was proposed by [16]. However, this method was recently 
looked at in great detail by [37], for the most optimistic 
case - dark matter in real space. There it was shown that, 
although the method accounts for broadband tilting of 
the underlying power spectrum, it was unable to correct 
for the shift of the peak due to mode-coupling. Since 
these mode coupling terms are even more enhanced in 
the halo- halo power spectrum (see [35]) we expect that 
this method will not correct for all of the physical shifts 
found in the previous section. 

In passing, we note that it is not straightforward to 
draw a direct connection between what we call the phys- 
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ical shift and what [37] describe as mode-coupling shifts. 
It is likely that what we have called a physical shift is an 
underestimate of the mode-coupling effects. To under- 
stand why, we note that whilst the Rcnormalized Per- 
turbation Theory formalism has yet to be extended to 
Haloes, we may suppose that the RPT decomposition of 
the 2-pt clustering signal will still be valid, i.e. there 
is some propagator which multiplies the linear theory 
power, and some sum of mode-coupling terms. Suppose 
the halo propagator has almost the same form as the 
dark matter propagator, so the effects of the non-linear 
bias mainly affect the mode coupling pieces. Since the 
propagator is akin to a Gaussian smoothing term, this 
term acts just like the simple Gaussian smoothing model 
we fit to the simulation data. If the effects of the other 
(mode-coupling) term were negligible, or did not depend 
strongly on halo mass, then we would expect the scale 
of the best fit smoothing filter to also be independent of 
halo mass. It is not, suggeting that the mode-coupling 
term depends on halo mass. That the scale of the best- fit 
smoothing filter is larger for the more massive halos (Ta- 
ble II and III) indicates that our best-fitting Gaussian 
filter is trying to account for some of the shifting that is 
actually due to the mode-coupling terms. 

We shall now pursue an analytic approach that we 
think provides insight in to the physics behind the shifts. 



IV. A PHYSICAL MODEL FOR THE SHIFTS 

This section presents a simple physical model for esti- 
mating the effects of nonlinear clustering and bias on the 
position of the local maximum of the correlation func- 
tion. [37] discuss a more accurate model for the corre- 
lation function of the dark matter; the approach below 
allows one to address how the peak shifts are affected if 
the measured correlation function comes from a biased 
tracer of the dark matter field. 



A. The pair conservation equation 



equation [58-61]: 
d[l + £(r,r)] 



v • 



[1+Z(t,t)]v 12 (t,t) -0, (21) 



where the divergence is with respect to the vector that 
separates the pair r = x\ — x 2l and the pairwisc infall 
velocity is 



«i2(r,r) 



((l + ft)(l + a 2 )(tti -tt 2 )) 
+ 



(22) 



where by statistical isotropy we used that (SiVi) = 
(S 2 v 2 ) = 0. 

We can rewrite Eq. (21) in a more convenient form, 
by changing time variable from conformal time r to the 
linear growth factor D + . In particular, if 



i] = lnD+, 



(23) 



then dr = dr]/(H,f), where H = dlna/dr and / = 
d In D + /d In a. We may also write velocities in a similar 
fashion and scale out their dependence on linear theory. 
Namely, v = —Hfu, where V • u = S in the linear theory. 
Then, dividing Eq. (21) by [1 + £(r,r)] yields, 



01n[l + g(r,q)] 
dr] 



«12 • V ln[l + £(r, 77)] = V • 1*12(1-, rj) . 

(24) 

Owing to the fact that large-scale flows have no vor- 
ticity, the pairwise velocities are directed along the sep- 
aration unit vector f , so U12 = u\ 2 r. Hence Eq. (24) 
becomes, 



d\n{l +gr,r 1 )} 
dr] 



"12 (r, n) 



dlnjl + gr,ri)}) 
dr 



e(r,r?) , 
(25) 



where we have defined 8(r, 77) = V • [1*12 (r, r?)f] to be the 
divergence of the pairwisc infall velocities ui 2 (r) . Note, 
that this equation may be thought of as a differential 
equation for ln(l + £) given an ansatz for U12 [59], or 
'vice- versa' [61]. 



The perturbed continuity equation for the collisionless 
CDM fluid can be written, 



a[i + a(x,T)] 

dr 



+ V 



(1 + <J(x,t))v(x,t) -0. (20) 



where5(x,r) = [p(x,r) — pb(r)]/ pb(r) , is the dimension- 
less density perturbation at comoving position x and con- 
formal time r (dr = dt/a(t), where a(r) is the expansion 
factor from the Friedmann equation); /^(t) is the homo- 
geneous background density; and u(x, r) = x' = dx/dr 
is the proper peculiar velocity field [57, 58]. 

We can now use Eq. (20) at position 1, say, multiply by 
(1 + S 2 ) for position 2, and add the same expression with 
indices 1 and 2 interchanged [6i = #(x,)]. Taking expec- 
tation values of the result yields the pair conservation 



B. Solution by characteristics 

The general solution of Eq. (25) can be found by the 
method of characteristics (see for example [62]), which 
illustrates quite clearly how any feature in the correlation 
function will move as clustering develops. 

The continuity equation (and thus the pair conserva- 
tion equation) is a prime example of a hyperbolic partial 
differential equation. Information propagates from the 
initial conditions to the final conditions through curves, 
called characteristics. The characteristics are simply the 
equations of motion of pairs, 



dr 
drj 



= -u 12 {r,r]). 



(26) 
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The solution of this equation gives r(r]), and this con- 
verts the left hand side of Eq. (25) into a total derivative. 
Thus, one obtains an ordinary differential equation along 
the characteristics: 



rflnjl + £(r,7 ? )] 
dr/ 



= 9(r,7 7 ) 



(27) 



and it should be understood that it is a function of time r\ 
only, after using the characteristic solution r(rj), Eq. (26). 
Thus Eq. (27) simply gives the logarithmic rate of change 
of the two-point correlation function as it evolves along 
the characteristic curve. The fully evolved correlation 
function may then be obtained straightforwardly, at any 
chosen epoch, through integration along the characteris- 
tic between the initial and final epoch: 



i + £M) = (i + £ohM)]) 



x exp 



e[v(r,77),7/]d7/ , (28) 



where rg (r, rj) is the initial separation that corresponds 
to r at time 77, and similarly for ry. The exponential 
factor comes from the fact that the correlation function 
is not conserved along characteristics because the right 
hand side of Eq. (27) is non-zero. Since we are mostly 
interested in significant growth after the initial pertur- 
bations are laid down (77 ^> 770), the term in the first 
parenthesis can be safely approximated as unity. Hence, 
all the evolution is encoded in 8 and the characteristics. 
Note that this solution is exact; it only becomes useful, 
though, if one can model the pairwise infall velocities. 
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FIG. 7: The flow of characteristics in linear theory for ini- 
tial separations close to the acoustic peak of the two point 
function, every 1 h~ x Mpc. For each scale we show results 
for dark matter (solid, solutions of Eq. 32), and linearly bi- 
ased tracers (solutions of Eq. 50) having z = bias factors of 
b = 1.4 (dashed) and b — 2 (dot-dashed). The peak in the 
linear correlation function is located at r = 106 h~ x Mpc for 
the cosmological model we use in this paper. 



C. Linear theory velocities 

For what follows, it will be convenient to define 

Uro,Vo) = e 2 " I J ^ Ji(kr)d"k, (29) 

where Po(k) is the power spectrum at some initial time 
770 = and where ji (y) = [sin(?/) — j/cos(y)]/y 2 is the first 
order spherical Bessel function. In linear theory, pairwise 
infall velocities, at time 77, can be written [58, 60] 

u n {r, n ) = 2 e 2 " J ^lj 1 (kr)d 3 k=^-e^Ur) ■ 

[63] . The divergence of pairwise velocities in linear theory 
can be obtained directly from Eq. (30) by taking the 
divergence, 



9(r,77) = V r • [ui 2 (r)f] 



l__d 

r 2 dr 



r 2 ui 2 {r)} 



2e 2ri J Po(k)j Q (kr)d 3 k , 

2e 2 "£o(r), (31) 



with £0 the initial (linear) correlation function at 770 = 0. 

Eq. (30) allows us to solve for the characteristics in 
linear theory. Two-point information at separation 7*0 
and time 770 = propagates by time 77 to a separation r 
(less than rg, due to clustering) so that, from Eq. (26) 



3 dr 



(32) 



Figure 7 shows the solution of this equation (r as a 
function of redshift) for initial separations r close to the 
acoustic peak of the two-point correlation function. If 
this were the only effect, i.e. if the right hand side of 
Eq. (27) were zero, then the correlation function would be 
conserved along the characteristics (solid blue line shown 
in Fig. 7) and this alone would give about 0.2% shift in 
the acoustic peak position by z = 0. However, as men- 
tioned above, the correlation function grows along the 
characteristics. This growth is governed by the diver- 
gence of the infall velocities, and, for large 77, it is this 
contribution which dominates. Indeed, we have not yet 
even included the linear amplification of the correlation 
function, resulting from the right hand side in Eq. (27). 

Including the divergence of infall velocities using 
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Eq. (31), makes Eq. (28) for the two-point function 

i + £M) = (i + £o[roM)]) 

x exp 



2 / eoMr.rOJe^dj/ 
o 



(33) 



If the flow of characteristics caused by the nonlinear term 
in Eq. (25) is ignored, then r « ~ ro, and so 



l + £(r) 



(l + &(r)) exp ^ (r)(e 2 "-l) 
l + &(r) c 2 ". 



(34) 



The final expression follows if the term in the exponen- 
tial is small; notice that it equals the linear perturbation 
theory expression for £ at time r\. 

At first sight, the solution of Eq. (33) appears to re- 
quire many evaluations of Eq. (32). However, the integral 
over rf in the exponential piece of Eq. (33) may be trans- 
formed using the characteristic curve, whence 



2e 2 " = d(e 2v - 1) = - 



3 dr 



(35) 



Thus on performing this change of variables, the term in 
the exponential of Eq. (33) becomes 



., r° dv 

r> f (r') 



(36) 



However, on noting that <i[r 3 £(r)] /r 3 = 3^(r)dr/r, we 
find that this may be further simplified to 

yrSfoCm) dx 

Jr 3 £ (r) 

Therefore, Eq. (33) is really rather simple: 



(37) 



1 + = (l + &(ro)) 



(38) 



and a single evaluation of Eq. (32) gives ro(r,rj), and 
hence the nonlinear value of £(r). 



D. Connections to previous work 

At late times e 2r ' 3> 1. Hence, on the large scales where 
£o < 1, Eq. (34) implies that 1 + £(r) w exp[S„(r)], 
where S^(r) = e 2l7 £o( r ) is the linearly evolved correla- 
tion function. This is precisely the relation between the 
correlation function of a lognormal field and that of the 
underlying Gaussian field from which it was derived. Of 
course, this analysis has assumed that r s» r v > « ro; 
Figure 7 shows that this is inappropriate at late times. 
Nevertheless, it provides a nice illustration of why the 
Lognormal has proved to be such a useful approxima- 
tion, and why the approximation breaks down [75] . Note 
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FIG. 8: The divergence of pairwise infall (dark matter) veloci- 
ties as a function of scale measured in numerical simulations 
at 2 = (solid squares) and at z — 1 (solid triangles). In lin- 
ear perturbation theory, 0/2 is equal to the linear correlation 
function - and at z — and z = 1 the linear models are rep- 
resented by solid red and blue lines, respectively. Notice that 
at z = 1 the acoustic peak is visible in O, but by z = the 
nonlinear effects have completely washed out any trace of it. 



that, both in linear theory and in the Lognormal approx- 
imation for the nonlinear evolution, the position of the 
acoustic peak does not shift [76] . 

Our Eq. (38) has the flavor of an approach pioneered 
by [59, 65], who argued that 



1 + £M) = {ro/rf 



(39) 



should provide a good approximation to nonlinear evolu- 
tion. In effect, their approach sets 



1 



ro\ 3 d In r 
r / dlnr 



(40) 



(41) 



•f(r,J7) 

If we set 1 + Co( r o) ~~ * F then we have 

1 + £(r,7i) = I — ) TTT ; 

our expression follows from inserting the linear velocities 
in the characteristics — it is not an ansatz. Note that this 
relation changes if nonlinear velocities are used. 



E. Inaccuracy of linear theory velocities 

In linear theory the divergence of infall velocities 
0(r, rf) is, modulo a factor of two, given by the linear 
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two-point function itself (c.f. Eq 31) and this has a 
static (independent of 77) peak at the unperturbed po- 
sition. Hence, there is a competition between 0, which 
prefers the peak to stay unshifted, and the flow of char- 
acteristics, which induce a shift towards smaller scales 
(Fig. 7). A consequence of this is that, using linear ve- 
locities is expected to underestimate the true peak shift. 
(The top left panel of Fig. 9 shows this explicitly, as we 
discuss later.) Using Eq. (38) one obtains a shift of about 
0.1% at z — 0, half of that due to the flow of character- 
istics. 

This underestimate results from the fact that, whilst 
the pairwise infall velocity may be reasonably well de- 
scribed by linear theory on large scales, its divergence 
deviates from linear theory more strongly, due to the 
scale dependence of nonlinear corrections [42, 66]. This 
is graphically illustrated in Fig. 8. Although 6 oc £ in 
linear theory, by z — 0, nonlinear effects have washed out 
any sign of an acoustic peak in 6! 

In practice, a characteristic that probes scales slightly 
smaller than the unperturbed acoustic peak will experi- 
ence more growth of the two-point function at late times. 
This leads directly to an enhancement that dominates 
over the effect of the flow of characteristics, and results 
in a substantially enhanced shift over the linear case (and 
as we will show this enhancement is about one order of 
magnitude). In this sense the flow of characteristics only 
gives a lower bound to the shift in the peak position due 
to mode coupling. Clearly, in order to proceed, we re- 
quire a model for the nonlinearity of the infall of pairwise 
velocities, and in particular its divergence G. 



F. Beyond linear theory velocities 

There are two types of nonlinear contributions to the 
pairwise infall velocity. This can be seen more clearly by 
rearranging Eq. (22) into the form, 



Eq. (42) thus leads to the following decomposition 



U12 



((Si + £ 2 )(«1 - U 2 )) + (<M 2 (ttl - u 2)) 
1+C 



(42) 



If we insert the standard perturbation theory (hereafter, 
PT) expansions for 5 and u [57], then we see that the 
first term in the numerator is second order in <5(x, 770), 
and the second term is of third order, which in linear 
theory averages to zero. We can set the denominator 
equal to unity, since £ is of order 10 -3 on the scales of 
interest and we are after much larger (1 — 10%) effects. 
As mentioned in the previous subsection, the effects from 
nonlinear mode-coupling on Ui 2 (r, 77) on these scales are 
negligible (<~ 1%), and hence play almost no role in shap- 
ing the characteristic curves (which, as we said, lead to 
shifts of only ~ 0.2% in linear theory for dark matter). 
They do, however, have a significant impact on the source 
term in the right-hand-side of Eq. (27), which dictates 
how fast the two-point function grows along the charac- 
teristics. 



0(r, n) = e 2 (r,r)) + Q 3 (r,r]), 



(43) 



where the two terms on the right-hand-side are defined 
02 = 2V • (S1U2) and ©3 = 2V • (S1S2U1). Considering 
the first term, on using the standard PT expansions for 
the density and divergence of the velocity field ([57] and 
see also footnote [63]), we find that ©2 can be written 

©2(r,r?) - 2 J P 59 (k,r ] )j Q (kr)d 3 k , (44) 



and 



P s %k, V ) = e 2 " P* e (k) + e 4 " < op (fc) , (45) 



is the cross-power spectrum of the density and veloc- 
ity divergence expanded to fourth order in the standard 
PT. The first term is the usual one from linear theory 
Pq 6 = Po, and Pfi oop is the 'one-loop' correction to P se 
from PT. The middle panel of Fig. 6 in [66] shows that 
this term describes rather well (much better than for the 
density power spectrum) the deviations from linear the- 
ory at large scales. Thus, 



Q 2 (r,n) = e 2 (r,n) + e 2 Ao °V(r,r 1 ). 



(46) 



Considering the second term in Eq. (43), we find that 



3 (r,r?) = 2V-<5i5 2 «i> 

x B m (k u k 2 ) , 



2e 4 " / d^Ae^ 1 - 

J M 



(47) 



where r = x\ — x 2 , k\ 2 = fei + k 2 and B ses 
is the density-velocity divergence-density bispectrum: 
(6(k 1 )6(k 2 )6(k 3 )) = B 5e5 (fe 1 ,fe2)5 D (fei +k 2 + k 3 ). Ap- 
pendix A provides explicit expressions for ©2 loop (r) and 
© 3 (r) expressed up to 1-Loop in the standard PT, and 
written in terms of the initial power spectrum. 

A substantially improved model for the nonlinear cor- 
relation function £ results from including these nonlinear 
terms in Eq. (28). Before showing this explicitly, the 
next subsection discusses how the effects of galaxy/halo 
biasing can be included in our analysis. 



G. Extension to biased tracers 

The analysis above has been useful for understanding 
the motion of the acoustic peak in the dark matter cor- 
relation function. However, since the observations will 
not measure the mass directly, but instead the cluster- 
ing of some set of biased tracers of the density field, i.e. 
some sampling of the galaxy distribution, the method of 
characteristic curve solutions will be more useful if we can 
extend it to describe these biased tracers. At first glance, 
it is not obvious that this can be done, owing to the fact 



17 



that halos, and the galaxies that they host, are created 
and destroyed through merging, so their comoving num- 
ber density is not conserved. Thus one might naively 
conclude that any such approach based on continuity ar- 
guments must be suspect. However, some thought shows 
that this problem is not insurmountable. 

Consider the motion of some halo today, its trajectory 
is the result of the previous history of motions of its con- 
stituent particles. Thus, for instance, one may speak of 
the motion of the center of mass of the particles that 
make up the halo, at, say, the present time. In partic- 
ular, one may also speak of the position and velocity of 
its center of mass even at high redshifts when the halo 
itself does not yet exist as a single virialized entity. This 
was the point made by [67]; provided appropriate care is 
taken of how the bias associated with these tracer parti- 
cles evolves, the continuity equation can indeed be used 
to relate £ to ui 2 . The argument above remains true if 
each halo is represented not by one but by many tracer 
particles, and the number of tracer particles depends on 
halo mass. The positions of each of these tracers can 
be followed back in time, so their number is conserved. 
These tracers have some effective bias factor at the time 
they are identified; provided one accounts for the evo- 
lution of this bias, the continuity equation can be used. 
Since the argument above works for any set of tracers, 
it is as valid for galaxies as for halos. Note in particular 
that detailed knowledge of the origin of the effective bias 
factor is unnecessary. E.g. if two sets of tracers have 
the same abundance and bias factor at one epoch, but 
one tracer populates a wide range of halo masses, and 
the other two narrow but rather separate mass bins, the 
evolution of the effective bias factor will be the same. 



Fortunately, describing the evolution of the bias for 
'objects' that arc neither created nor destroyed is rather 
straightforward [67-72] : For a set of tracer particles that 
are related to the underlying dark matter through a lin- 
ear, local, deterministic mapping, the time evolution of 
their bias (b(i]) = 6 a (x 7 r])/S(x,ri) where a represents ei- 
ther haloes or galaxies), can be written 

b( V ) -l = (h- l)e-" , (48) 

where hi denotes the bias at the initial time rj = 0. Thus 
to incorporate this bias model into our theoretical model, 
we must simply make the following replacements: 

& - &? Co ; e 2 -> b( v ) e 2 ; e 3 - b( v ) 2 e 3 , (49) 

in the expressions above. Here we have used the stan- 
dard assumption that the velocity field of any set of 
biased tracers is itself unbiased, and that O3 depends 
only quadratically on the density field, where we have 
neglected sub-leading terms (see [70]). 

With these changes, Eq. (32) for the characteristics 
becomes 

e 2 "-l + 2(6 i -l)(e"-l)= / _ . (50) 

Fig. 7 shows solutions to this expression for tracers that 
have bias factors of b = 1.4 and b = 2 at z = 0. It shows 
that the flow of characteristics towards small scales is 
enhanced if b > 1; and this is as expected, because infall 
velocities are proportional to the bias factor [67]. 

Our model for the nonlinear correlation function of bi- 
ased tracers means that Eq. (28) becomes 



£M) = (l + 6?&[ro(r,f?)]) 



x cxp 



drj'(b(r)') 6 2 [ry(r, 77), r/] + b{-q'f 3 [ry(r, f?),?/]) 



(51) 



Note that the linear theory solution of this equation may 
be recovered directly by setting: 2 = © 2 ; B3 = 0; and 
= r = r in the expression above. Whence, 

£(r)«&fa) 2 6,(r)e 2 " , (52) 
and this is the generalization of Eq. (34). 



H. Comparison with simulations 

Figure 9 compares our model for the nonlinear correla- 
tion function, Eq. (51), with our measurements of (real- 
space) £ for the dark matter (top) and halos (bottom) at 
z = (left) and z = 0.5 (right). The halo measurements 
are the same as those presented previously, except that 



now we only show scales which are within <~ 15 h 1 Mpc 
of the initial acoustic peak. 

Our model for the dark matter, Eq. (28), matches 
the measurements rather well; the solid lines are a sub- 
stantial improvement over linear theory (dashed). Our 
model predicts that the peak has shifted to 105 Mpc 
by z — 0.5, about a one percent effect; this is in good 
agreement with a more rigorous calculation based on 
RPT [37]. By z = our model for predicts that 
the peak has shifted to 98 ft -1 Mpc, roughly an ~ 8 
percent shift! This disagrees by over a factor of four 
with the RPT calculation (apparent and physical shifts 
included); this overshoot is not surprising, when given 
the fact that one-loop PT is known to overestimate the 
nonlinear power spectrum by tens of percent on small 
scales, even though the one-loop density-velocity diver- 
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FIG. 9: The real-space two-point correlation function for dark matter (top) and halos (bottom) at z = (left) and z = 0.5 
(right). Table I describes the three bins in halo mass; £hh is largest for the most massive halos. The dashed lines in the top 
panels show linear theory for the dark matter, solid lines are the predictions of our model, Eq. (51), and symbols show the 
measurements. Dotted line in the top left panel shows our model when only linear theory velocities are used; it is almost 
indistinguishable from simple linear theory, demonstrating that inclusion of the nonlinear contributions to the (divergence of 
the) velocity field is vital. 
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gence power spectrum does well at reproducing the cross- 
power spectrum as measured from numerical simulations 
at intermediate scales [66]. 

Turning now to the results for the dark matter halos, 
we see that Eq. (51) provides a very good description of 
the measurements. We emphasize that there are no free 
parameters in this model. The only non-cosmological pa- 
rameters in the model are the bias factors and as dis- 
cussed earlier, these are measured directly from the sim- 
ulations to make the predictions (see Section III B for our 
estimated values for the halo mass bins listed in Table I). 

When the bias factor is large, then the dominant non- 
linear correction comes from 63 because it scales as b 2 . 
For the dark matter, the nonlinear correction coming 
from 02 loop is the dominant one. The figure shows that 
our model does not predict any significant trend with 
halo mass, although this would likely change if we were 
to include nonlinear bias (e.g. [35] suggest higher mass 
halos will show enhanced nonlinear effects). Our model 
requires knowledge of how these nonlinear bias terms 
evolve (i.e., the analog of Eq. 48): this evolution is given 
in [73, 74]. 



V. CONCLUSIONS 

We have used analytic methods and a very large ensem- 
ble of numerical simulations to study how the position 
of the baryonic acoustic peak in the two-point correla- 
tion function, £, remnant of the tight coupling between 
photons and baryons before recombination, is affected 
by the clustering systematics: nonlinear mass evolution, 
bias and redshift space distortions; and we have examined 
these effects as a function of cosmological epoch and as a 
function of several trace particle types - i.e. halo samples 
picked to evolve with constant comoving number density. 

We have investigated a toy-model for the evolution 
of £ (Sec II) that was simply a Gaussian bump plus a 
power-law and we showed that, if nonlinear evolution 
was manifest as a Gaussian smoothing of the true £, then 
the acoustic scale was not well recovered through simply 
measuring the local maximum - and this we described 
as an apparent shift of the peak. However, if there was a 
change in the underlying broad band power then this we 
said would lead to a a physical motion of the peak. 

We presented results from our numerical simulations 
(Sec III A). Our total simulated volume corresponded to 
~ 105 Gpc 3 /i -3 , approximately the same size volume 
that the proposed Stage IV, JDEM mission, ADEPT in- 
tends to survey [43]. Therefore our results and analysis 
arc of direct relevance to that and similar missions. From 
these simulations we measured £ for the dark matter and 
haloes. We found, at z = in both real and redshift 
space, that the true position and shape of the linear the- 
ory function did not match well that of the measured 
data - there being an enhanced signal on scales smaller 
than the unperturbed peak scale. 

We then performed a more careful analysis, and fit- 



ted the correlation function data using the Gaussian 
smoothed linear theory model. This provided a much 
improved fit. In all cases the inferred peak positions 
from these models were shifted to smaller scales, with 
typical shifts being of the order ~ 0.5 — 3.0 Mpc 

including a factor of ~ 2 correction for the transfer 
function; the shifts were enhanced for the the highest 
mass haloes/rarest objects and in redshift space. How- 
ever they were alleviated for higher redshifts. We con- 
cluded that this was direct evidence for 'apparent motion' 
of the acoustic peak. We also noted that many of the re- 
cently proposed BAO reconstruction methods do attempt 
to account for this apparent shifting of the peak. 

We then showed that the smoothed linear model was 
not a perfect fit to the data, in particular for highly bi- 
ased haloes and their galaxies the fit was poor. Using 
the x 2 test we ruled out this model at the 99% signifi- 
cance. Furthermore, through inspection of the residuals 
of the fitting we found ~ 10 — 20% excess of amplitude 
on scales smaller than the unperturbed acoustic scale. 
We concluded that this was supporting evidence for the 
hypothesis that non-linear evolution was inducing a phys- 
ical motion of the acoustic peak. We characterized the 
physical shifts by finding the local maximum of smooth 
polynomial fit to the data and subtracting from it the lin- 
ear model peak position. We found that these shifts were 
of the order ~ 0.0 — 3.0 h^ 1 Mpc for the samples consid- 
ered. We noted that these - which represent broad band 
tilting plus the more pernicious mode coupling effects - 
are not accounted for in recently proposed schemes to 
correct the signal in the power-spectrum. 

In our analysis of the simulation data we also presented 
evidence that the simple Gaussian-based calculation for 
the variance (Eq. 13) of £ that included the Poisson shot- 
noise contribution provided a good description of the ex- 
pected error on the measured £ for haloes (Figs. 4 and 5). 
In detail the Gaussian error model was found to under- 
estimate the simulations for haloes at the present day. 
Adding non-Gaussian terms from the trispectrum and 
bispectrum may improve this further. 

For future BAO missions that aim to use the power 
spectrum of clusters, the expected sample variance esti- 
mates that use the Gaussian plus Poisson model, will give 
reasonable estimates of the variance. Directly extrapo- 
lating our analysis to make a statement about the vari- 
ance estimates for BAO galaxy surveys is complicated. 
Our analysis has dealt with the clustering of the halo 
centers and does not include the virial motions of parti- 
cles/galaxies internal to the halo - adding in this layer 
of reality may lead to increased variance. Furthermore, 
if galaxies appear only in haloes, then using the galaxy 
number density estimate as the Poisson shot noise er- 
ror as is common practice, will underestimate the sample 
variance when there is more than one galaxy per halo. 
We expect that the effective number density of the haloes 
that host the galaxies in the survey will be a better reflec- 
tion of the errors. We shall reserve this issue for future 
study. 



20 



We then presented an analytic model that was able 
to capture the main observed affects from the non-linear 
evolution of the mass and bias. The model was based 
upon a study of the gravitationally driven mean stream- 
ing motions of particle pairs. These motions both smooth 
out the initial peak, and, more importantly, shift it 
(Fig. 7). In essence, our model simultaneously accounts 
for both the smoothing and the shifting of the acoustic 
peak. We first discussed the model in the context of the 
dark matter (Eq. 28), and then showed how it could be 
extended to describe the nonlinear evolution of £ for bi- 
ased tracers, such as galaxies and clusters of galaxies as 
well (Eq. 51). For the dark matter, our approach is less 
reliable than that of RPT (see [37] for a discussion of 
this). However, we think it has substantial merit, owing 
to the fact that it permits a simple description of how the 
shifting of the acoustic peak is modified for biased tracer 
particles. It also allows us to see the problem from a 
different perspective. One could combine the strength of 
both methods, by replacing the modeling of the diver- 
gence of pairwise velocities by its RPT description, for 
that one would need to calculate the bispectrum contri- 
bution to 03. 

The measured shifts in the acoustic peak position for 
the dark matter and the biased tracers, are qualitatively 
consistent with the effects of the clustering systematics on 
the power spectrum [35]. This owes to the fact that the 
power spectrum and correlation functions are a Fourier 
transform pair. Thus small scale damping and tilting of 
the linear power spectrum can lead to both smoothing 
and tilting of the correlation function, and the genera- 
tion of the measured apparent and physical motion of 
the peak. 



However, these recovered shift values appear substan- 
tially larger than those currently quoted in the literature 
from other analytic arguments [30]. One possible expla- 
nation for this is that the divergence of the pairwise veloc- 
ity field is substantially more non-linear than the density 
field on these large scales (Fig. 8). Had we simply used 
linear theory velocities in our analytic model then we 
would have considerably underestimated the measured 
shifts. Using perturbation theory was crucial (Eqs. 43- 
47) for our model to get the close agreement with the 
numerical measurements. 

If unaccounted for, the percent level changes we have 
measured in the acoustic scale will lead to biased de- 
terminations of cosmological parameters (and see [11]). 
However, the agreement between our model and the sim- 
ulations suggests that, although such pernicious shifts 
are present, it may be possible to construct analytic tools 
that allow us to correct for them. This is the subject of 
ongoing work. 
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APPENDIX A: THE DIVERGENCE OF INFALL VELOCITIES IN PERTURBATION THEORY 

This Appendix provides expressions for 92 loop (r) and 3 (r) from the standard PT. 



1. e" 0015 in the standard PT 



O;!, loop is given by Eqs. (44) and (45), and is an integral over the 1-Loop contribution to the velocity divergence- 
density power spectrum. In the standard Perturbation Theory this can be written [66]: 



AioopW = 2 / F 2 (k - q, q) G 2 (k - q, q) P (\k - q\) P (q) <Fq + 3 P (k) [F 3 (k,q) + G 3 (k,q)} P (q) d 6 q, (Al) 



where the functions F 2 (k,q) and G 2 (k,q) are the second order, symmetric, density and velocity divergence kernels 
from PT [57]. These are written: 



5 1 



F 2 (k,q) = 7 + 2^(- + f)+ 7 , 



^ /. n 3 1 (K q\ 4 2 



k , <T 

q 



(A2) 
(A3) 
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where Hk,q = k ■ q/\k\\q\. The functions F 3 (k,q) and G 3 (k,q) are the angle averages of the third order PT density 
and velocity kernels. These may be written: 



6fc 6 - 79fcV + 50k 2 q 4 - 21q 6 (q 2 - k 2 f(7q 2 + 2k 2 ) , 
1 „ ; In 



63fcV ' 42fc 3 <? 5 

6fc 6 - 41fcV + 2fcV - 3 g 6 (q 2 - k 2 ) 3 (q 2 + 2k 2 ) 



k + q 



21/cV 



+ 



Uk 3 q 5 



In 



k — q 
k + q' 



k — q 



;(A4) 
(A5) 



2. 6 3 (r) in the standard PT 

83 (r) is related to the density- velocity divergence-density bispectrum through two Fourier transforms (Eq. 47). In 
the standard PT this bispectrum is: 



B 5e5 (k u k 2 , fe 3 ) = 2F 2 (k 2 , k 3 )P (k 2 )P (h) + 2G 2 (k 1 ,k 3 )P {k 1 )P (k 3 ) + 2F 2 {k u k 2 )P (k 1 )P (k 2 ) 



(A6) 



In order to proceed we require some further pieces of information. Firstly, the closure relation for fc-modes gives 
us k 3 = — fei — k 2 . Secondly, statistical homogeneity and isotropy means that the bispectrum can be written as a 
function of three variables: the length of two sides of a triangle and the angles between them, i.e. we should at the 
end of our calculation be able to write B 5es (ki,k 2 ,k 3 ) = B S0S (ki,k 2 , #12). Thirdly, the addition theorem for spherical 
harmonics allows us to re-write the angles between any two vectors in terms of their own angles in some arbitrary 
Cartesian system: 



cos 0i2 — cos 0i cos 2 + sin 0\ sin # 2 cos(0i — <fo) 



(A7) 



where the angle between the two vectors k\ \k\,0\,§\) and fc 2 {k2,02,<p2} is #12- Some lengthy algebra then leads 
us to the following expression for 83 (r): 



8 3 (r) 



d d kP Q 



«*)^) fc /^)[(!*-* +1 )i^r-5(? + F 



*§(r)*rV) + *2 2 (r)*{(r)\ + - [^(r^rV) + *o 
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-* (r)*l\r) 
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^(r)*rV) 



where we have introduced the useful auxiliary function 

*rW= / d 3 qP (q)je(qr) q" 



(A8) 



(A9) 
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